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Accurate modeling of left-handed media using finite-difference time-domain method 
and finite-size effects of a left-handed medium slab on the image quality revisited 
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The letter contains an important message regarding the numerical modeling of left-handed media 
(LHM) using the finite- difference time-domain (FDTD) method which remains at the moment one 
of the main techniques used in studies of these exotic materials. It is shown that conventional 
implementation of the dispersive FDTD method leads to inaccurate description of evanescent waves 
in the LHM. This defect can be corrected using the spatial averaging at the interfaces. However, a 
number of results obtained using conventional FDTD method has to be reconsidered. For instance, 
the accurate simulation of sub-wavelength imaging by the finite-sized slabs of left-handed media 
does not reveal the cavity effect reported in [Phys. Rev. Lett. 92, 107404 (2004)]. Hence the finite 
transverse dimension of LHM slabs does not have significant effects on the sub-wavelength image 
quality, in contrary to previous assertions. 
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The finite-difference time-domain method (FDTD) is 
known as one of the powerful numerical techniques in 
electrodynamics [l|. Being simple in implementation it 
has been proved to be very popular among researchers. 
This method is assumed to be extremely accurate since it 
involves direct numerical solution of Maxwell equations 
which are known as the basis of classical electrodynamics. 
However, the implicit belief in FDTD sometimes results 
in attributing certain physical properties to some elec- 
tromagnetic structures based on simulation results. One 
typical example is a long row of works on FDTD model- 
ing of the left-handed media (LHM), materials with neg- 
ative permittivity and permeability 2]. Such materials 
are not yet available experimentally and thus numerical 
simulations still remain one of the most common ways 
in exploring their properties and applications. However, 
as it will be shown in this letter, the conventional imple- 
mentation of FDTD for modeling of LHM in the same 
manner as for usual dispersive dielectric materials leads 
to incorrect simulation results. It may seem that the 
conventional FDTD has been verified in the literature: 
the negative refraction effect which is inherent to the 
boundary between the free space and LHM was observed 
and the planar superlens behaviour has been successfully 
demonstrated 0, fl • Actually, this only means that 
the LHM is correctly modeled for the case of propagat- 
ing waves. As soon as the evanescent waves are consid- 
ered the conventional implementation of the FDTD fails. 
Usually, the evanescent waves decay exponentially over 
the distance and thus they are concentrated in the close 
vicinity of sources, that is why conventional FDTD mod- 
eling of usual materials does not suffer from this trouble. 
In the case of LHM, the evanescent waves play key roles 
and have to be modeled accurately because of the perfect 
lens effect 6] . A slab of LHM effectively amplifies evanes- 
cent waves which normally decay in usual materials and 
allows transmission of sub-wavelength details of sources 
to significant distances. Ideally lossless LHM slabs pro- 
vide unlimited sub-wavelength resolution. However in 



realistic situations, the resolution is restricted by losses 
and the thickness of the slab [jj , as well as the mismatch 
between the LHM and its surrounding medium 

In order to illustrate what we mean by incorrect de- 
scription of evanescent waves provided by conventional 
implementation of FDTD, we have simulated propaga- 
tion of plane electromagnetic waves (with TM polariza- 
tion and various transverse wave vectors) through an in- 
finite slab of LHM in the free space. We choose this 
problem since an analytical expression for transmission 
coefficient through the LHM slab is available p| , that al- 
lows us to check validity of our numerical results. The 
LHM slabs with relative permittivity and permeability 
£ = [i = —I — O.OOlj and thickness d = A/5, where 
A is wavelength in the free space, are tested. The fre- 
quency dispersion of the LHM has been modeled by the 
Drude model with the plasma frequency u p = \J2uj and 
the collision frequency 7 = 0.0005cj, where uo is the op- 
erating frequency. It is implemented in FDTD using the 
auxiliary differential equation (ADE) method 1]. The 
two-dimensional (2-D) simulation domain is bounded by 
perfectly matched layers (PMLs) and periodical bound- 
ary conditions (PBCs) as illustrated in Fig.^a). A soft 
current sheet source (which allows scattered waves to 
pass through) with phase delay corresponding to differ- 
ent transverse wave vectors k x placed at the distance d/2 
from the slab is used as excitation. The PBCs are also 
specified by the particular wave vector k x . The source 
is slowly and smoothly switched to its maximum value 
in order to avoid exciting other frequency components 
The computation continues until the steady state is 
reached. The simulation has been done for both propa- 
gating (k x < /c, where k is wave vector in the free space) 
and evanescent (k x > k) waves. The Berenger's original 
PMLs 9] are used for absorbing propagating waves, and 
the modified PMLs ^} are applied when the transmis- 
sion coefficient for evanescent waves is calculated. The 
PMLs are placed at A/2 distance away from the slab. We 
use Yee's square grid with periods Ax = Ay — A/100 
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FIG. 1: (Color online) (a) The simulation domain for calcu- 
lation of the transmission coefficient, (b) The comparison of 
transmission coefficients as functions of transverse wave vec- 
tor k x for the LHM slab with e — fi — — 1 — 0.00T/ and thick- 
ness d = A/5, calculated using conventional FDTD method 
(circles) with the analytical result (solid line). The result pro- 
vided by the scheme with averaging at the boundaries is given 
by crosses. 

and the time step At = Ax/V^c, where c is the speed of 
light in the free space, chosen according to the Courant 
stability condition 1]. Since infinite structures can be 
truncated with any period, in order to save computation 
time, we use only four FDTD cells along x-direction. 

The calculated transmission coefficient from the source 
plane to the image plane (located at the distance d/2 
from the other side of the slab) as function of k x is pre- 
sented in Fig. H^b) by circles. The reference curve cal- 
culated using the analytical expression |g is also shown 
for comparison. It is clearly visible that the numerical 
results from the conventional FDTD are correct only for 
k x < 2ko. This range of k x covers all propagating waves 
(k x < ko) and a small part of weakly decaying evanescent 
waves (ko < k x < 2ko). For the evanescent waves with 
k x > 2ko the numerical results dramatically differ from 
the analytical results and the former shows resonant be- 
havior with a strong peak at k x = 2.4/co- This effect can 
be explained as resonant excitation of a 'numerical sur- 
face plasmon' at the back interface of LHM slab. Similar 
phenomena can be observed in the case of metallic slabs 
p or for unmatched LHM 8] , but in this particular case 
it is purely numerical artefact. The incorrect behavior 
of numerical solutions remains similar if the FDTD grid 
period is reduced to A/200 and A/400, but the resonance 
shifts to k x = 2.8/co an d k x = 3.2/co, respectively. 

The presence of 'numerical surface plasmons' provides 
evidence that the boundaries between the LHM and the 
free space have not been modeled accurately. If at the 
boundaries the mean value of permittivity of LHM and 
the free space is used for updating the tangential com- 
ponent of electric field (which is equivalent to the spatial 
averaging suggested in 11]) then the spurious 'numerical 
surface plasmons' disappear and the modeling happens 
to be extremely accurate. The transmission coefficient 
calculated using the proposed spatial averaging at the 
boundaries are presented in Fig. H^b) by crosses. It is 



clear that the result repeats the estimated analytical val- 
ues with very good accuracy for the whole spatial spec- 
trum of waves. The calculation has been performed for 
Ax = Ay = A/ 100, and it remains accurate even for 
larger grid periods. The above simple test allows us to 
conclude, that the conventional FDTD method fails to 
describe the propagation of high-order evanescent waves 
in the LHM if no corrections at the boundaries of the 
LHM are made. As a result, a number of previously ob- 
tained results using the conventional FDTD have to be 
reconsidered. This especially concerns the modeling of 
sub- wavelength imaging by LHM slabs Q which involves 
operation with evanescent waves. Note, that the simu- 
lations where only propagating waves are involved (e.g. 
demonstration of negative refraction for an obliquely in- 
cident plane wave) are not affected by this problem. 

The numerical transmission coefficient for LHM slabs 
reported in [l2| is a typical example of using the FDTD 
method without averaging. The study can actually be 
treated as an investigation of the 'numerical surface plas- 
mons', their sensitivity to losses and efficiency of their 
excitation for various thicknesses of the slab. Unfortu- 
nately, these results have no relations with the properties 
of actual LHM slabs. The performance of the LHM slab 
as a sub-wavelength imaging device indeed depends on 
the losses and the thickness of the slab as it is shown in 
[3 , but it is completely different dependence as compared 
to results reported in [l^ . 

One of the recent most puzzling results related to the 
quality of imaging provided by LHM slabs is reported in 
p^ . It is claimed that the operation of the finite-sized 
structures is significantly affected by their transverse di- 
mensions. Having the FDTD code with spatial averaging 
at the interfaces which has been proven to be accurate, 
we decide to check this statement. The finite-sized slabs 
of LHM excited by magnetic current sources are modeled 
for three different transverse dimensions: L = A, 2 A and 
4A, as illustrated by the sketch in Fig. The parame- 
ters of the LHM slab (e = fi = -1 - O.OOlj, d = A/5) 
and distance between the source and the front interface 
equal to d/2 are kept the same for all simulations. The 
simulation domain is truncated by PMLs located at A/2 
distance away from the LHM slab and both source and 
image planes. The same grid periods Ax — Ay = A/100 
and the time step At = Ax / \^2c are used as for previous 
plane- wave simulations. The source is switched slowly 
and smoothly in order to avoid contributions from unde- 
sired frequency components p| and the simulations last 
until the steady state is reached. The intensity distribu- 
tions in the image planes for all three cases of different 
transverse dimensions are plotted in Fig. It is clear 
that the image quality is practically unaffected by the 
transverse size of the slab. The distributions repeat the 
source distribution, which is plotted in the same figure, 
with good sub-wavelength resolution. We do not observe 
any distortion of images caused by the fmiteness of the 
transmission device, in contrast to the conclusions made 
in The slight disagreement between the image and 
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FIG. 2: (Color online) The magnetic field intensities at the 
image planes of the planar LHM lenses {e = \i — — 1 — O.OOlj, 
d = A/5) with various transverse dimensions: L = A, 2A and 
4A. 

the source is due to the finite resolution of lenses caused 
by the losses in LHM. We suppose that the resonant ef- 
fects and image distortions related to transverse dimen- 
sions reported in 13] can be interpreted as excitation of 
'numerical surface plasmons' at the interfaces of the slab 
and can be observed only in inaccurate FDTD simula- 
tions without spatial averaging at the boundaries. Thus, 
these effects are purely numerical and have no relation to 
the properties of actual LHM slabs. In reality, the imag- 
ing performance of finite-sized LHM slabs is unaffected 
by their transverse dimensions. This statement also has 
been confirmed by full-wave electromagnetic simulation 
using Ansoft HFSS™ package. 

In order to illustrate this statement we have performed 
the simulation for a LHM slab with transverse size L = A 
excited by two magnetic line sources placed at A/8 dis- 
tance between each other using the FDTD method with 
spatial averaging at the boundaries. The rest of param- 
eters is the same as before. The distance between the 
sources is larger than the resolution of the lens which 
should be better than A/ 12 based on the fact that the 
transfer function plotted in Fig. ^ is close to unity for 
k x /k < 6. This allows us to expect two well-resolved 
maxima in the image plane. The distribution of mag- 
netic field intensity in a sub-domain near the LHM slab 
(the actual FDTD domain is larger) in the steady state 
is presented in Fig. Efa). The distribution in the image 
plane is shown in Fig.^c) together with the source. Two 
maxima at the distance of A/8 are clearly visible in the 
image plane. This confirms sub- wavelength imaging ca- 
pability of the LHM lens with only one wavelength width. 
Thus, we do not observe any limitations on the function- 
ality of the LHM slabs as sub-wavelength imaging de- 
vices due to their finite transverse dimensions. However, 
if the same system is modeled by FDTD method without 
corrections at the boundaries then completely different 
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FIG. 3: (Color online) The magnetic field intensity distribu- 
tions around the LHM slab with transverse size L — A excited 
by two magnetic line sources placed at A/8 distance between 
each other obtained using FDTD method (a) with and (b) 
without spatial averaging, (c) The image and source plane 
cuts. 



distribution of the field around the slab is observed, see 
Fig. [3b). The distribution of the field along the slab 
interface is smooth which once again confirms that high- 
order evanescent waves are not correctly modeled. As a 
result, the sub- wavelength details of the source are not 
resolved in the image plane. Instead of expected two 
closely located maxima the intensity distribution in the 
image plane has only one wide maximum, see Fig.|H{c). 

We suppose that the presented comparison between 
FDTD models with and without spatial averaging at 
the boundaries clearly demonstrates the limitation of the 
conventional FDTD method for modeling of LHM lenses. 
The significant discrepancies appear only in the cases 
when evanescent waves are involved. However, we en- 
courage to use the model with corrected updating equa- 
tions at the boundaries (with spatial averaging of per- 
mittivity as we propose or with averaging of the current 
as proposed in [ll[ which are equivalent) in all cases in 
order to avoid numerical artefacts. The spatial averag- 
ing is known as a second-order correction in the case of 
boundaries between dielectrics with positive permittiv- 
ity [l4j, but in the case of evanescent waves in LHM it 
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transforms into essential and mandatory correction. 

In addition to the corrections at the boundaries we 
would like to stress a few other aspects which are impor- 
tant to accurate FDTD modeling of LHM. The numerical 
dispersion is usually assumed to have very little influence 
on the quality of FDTD simulation. However, the proper- 
ties of LHM are known to be extremely sensitive to minor 
changes of material parameters. A typical example is the 
degradation of resolution due to small variations of ma- 
terial parameters from e = fi = — 1 0, 0- This means 
that even a tiny amount of numerical dispersion may 
cause crucial discrepancy between simulated and theo- 
retical results. Fortunately, the numerical dispersion can 
be easily analyzed analytically and an estimation for the 
difference between effective numerical material parame- 
ters and the 'real' ones are known This allows to ad- 
just parameters of the FDTD simulation to ensure that 
the effective numerical material parameters correspond 
to required values. We suggest using this adjustment in 
the cases of large time steps and small losses. 

It is well known that the switching time consider- 
ably influences the oscillation of images created by LHM 
lenses. The switching time equal to thirty periods was 
used in [1,0- However, perhaps this is the reason why 
no stable images could be obtained in 0, . Recently, 
it was reported in that a switching time equal to one 
hundred periods is required to obtain stable images. We 
have used such a switching time in all our simulations and 
we recommend to pay attention to the issue of switching 
time in all FDTD simulations of LHMs. The high-order 
evanescent waves travel very slowly in the LHM slabs and 
the procedure of the amplification of evanescent waves 
requires very long time to reach the steady-state. That 
is why, in addition to smooth and slow switching of the 
source we recommend to add losses into the simulations 
in order to limit the spatial spectrum of evanescent waves 
which are involved in operation. Otherwise, in the loss- 
less case the transient processes may last extremely long. 

The diffraction on the wedges and corners also may 
provide certain problems because of the singularity ef- 
fects reported in However, in the case of LHM with 
e = n « — 1 the singular behavior disappears and the 
wedges operate mainly as retroreflectors 17]. That is 
why, in the simulations presented in this letter we have 
not observed any singularities at the wedges. 

In conclusion, we have demonstrated in the letter that 
the conventional FDTD method for modeling of LHM 
leads to inaccurate description of high-order evanescent 
waves and does not simulate sub- wavelength imaging cor- 
rectly. The simulations suffer from the artefact of 'nu- 
merical surface plasmons' which appear at the interfaces 
of LHM. In order to solve this problem and ensure accu- 
rate FDTD modeling, a spatial averaging scheme at the 
boundaries has been proposed in the same manner as in 
[ll| . This technique has been tested on the analytically 
solvable example of a plane wave propagation through an 
infinite LHM slab and extremely good accuracy of simu- 



lation has been revealed for all angles of incidence includ- 
ing high-order evanescent waves. The finite-sized slabs of 
LHM excited by a single line source have been modeled 
using the proposed technique for various transverse di- 
mensions of the slab and the sub- wavelength imaging ca- 
pability of the structures has been confirmed. No distor- 
tions of the image due to the finite transverse dimensions 
of the structure have been revealed, in contrast to the 
results reported in . We suspect that the resonant ef- 
fects due to finite transverse dimensions of the LHM slab 
reported in ^| are caused by excitation of the 'numerical 
surface plasmons' and are completely numerical. These 
effects are absent in real structures and there are no re- 
strictions on functionality of LHM sub- wavelength lenses 
due to their transverse dimensions. Our FDTD simula- 
tions with spatial averaging at the boundaries demon- 
strate the imaging of two line sources located at sub- 
wavelength distance between each other by a LHM slab 
with only one-wavelength width. Keeping in mind that 
simulations using the FDTD method remain one of the 
most popular approaches in the studies of LHM, the gen- 
eral recommendations on how to accurately model LHM 
using the FDTD method are also provided. 
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